Measuring orbital interaction using quantum 

information theory 

Jorg Rissler^'*, Reinhard M. Noack^, Steven R. White ^ 

^Fachbereich Physik, Philipps-Universitdt Marburg, D-35032 Marburg, Germany 

^Department of Physics and Astronomy, University of California, Irvine CA 

92697-4 575, USA 

Abstract 

Quantum information theory gives rise to a straightforward definition of the inter- 
action of electrons Ip^q in two orbitals p, q for a given many-body wave function. 
A convenient way to calculate the von Neumann entropies needed is presented in 
this work, and the orbital interaction Ip^g is successfully tested for different types 
of chemical bonds. As an example of an application of Ip^g beyond the interpreta- 
tion of wave functions, Ip^q is then used to investigate the ordering problem in the 
density- matrix renormalization group. 
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1 Introduction 



Orbital interaction is used mainly as a qualitative concept in chemistry. Exam- 
ples are frontier orbital theory or the isolobal principle. Two orbitals are said 
to interact when they have similar energy and/or matching spatial distribu- 
tion and/or matching occupation. This concept is also important in theoreti- 
cal chemistry: for example, in complete-active-space calculations, one chooses 
(interacting) orbitals to form the active space. However, there is no unique 
quantitative measure of orbital interaction. 

One way to obtain a quantitative measure of orbital interaction is to utilise 
the concept of entanglement. Namely, if one divides the Hilbert space for a 
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given correlated wave function into two parts, then one can determine the 
entanglement of the two parts for the given wave function. For example, a 
division of a Hilbert space into one orbital and the space spanned by all the 
other orbitals defines the entanglement of this specific orbital with the rest 
of the Hilbert space. The more an orbital is entangled, the more it exchanges 
information with the other orbitals within the wave function. It is therefore 
natural to identify entanglement with orbital interaction. The advantage of 
this interpretation is that quantum information theory [1,2] defines a quanti- 
tative measure of the entanglement and hence of orbital interaction: the von 
Neumann entropy S. The disadvantage of using such a definition in calcu- 
lations is that is a function of density-matrix eigenvalues. Consequently, a 
correlated wave function must be obtained before a von Neumann entropy can 
be calculated. 

In order to proceed, one must use a method that produces an approximate 
wave function from which the von Neumann entropies can be calculated in 
a simple fashion. In addition, the approximate results should be close to 
those for the exact wave function. One suitable technique is to solve the full- 
configuration-interaction (FCI) problem using the density-matrix renormaliza- 
tion group (DMRG) [3]. This method produces wave functions and energies 
with a well-defined error. In addition to control of the accuracy of the wave 
function, the DMRG allows easy access to operators and expectation values 
in second-quantised form. We use this approach in order to calculate the von 
Neumann entropies in a straightforward way. 

In the following, we argue that a reliable measure of the orbital interaction 
between two orbitals can be formed by subtracting the von Neumann entropy 
of the two orbitals taken together from the sum of the individual von Neumann 
entropies of each orbital. We then explore the utility of this measure and how 
their properties depend on the accuracy of the wave function, we investigate 
four example molecules which contain different types of chemical bonds: LiF, 
CO, N2, and F2. It turns out that the results do not depend significantly on 
the accuracy of the wave function. It is therefore sufficient to determine the 
orbital interaction using a comparatively inexpensive calculation. 

We then apply this definition of orbital interaction to a problem in theoretical 
chemistry: namely, the ordering problem in the DMRG algorithm. One ingre- 
dient needed to define the DMRG algorithm is to order the orbitals that span 
the Hilbert space onto a one-dimensional lattice. This ordering has a signifi- 
cant infiuence on the convergence of the energy. In order to find an ordering 
that ensures a good convergence, the literature offers several approaches [4,5], 
which all have one aspect in common: they define a measure of orbital inter- 
action and order the orbitals in a way so that strongly interacting orbitals are 
near each other. Therefore, the measure of orbital interaction defined in this 
work opens up a new way of approaching this problem. 
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These topics are treated in the foUowing order: First, an introdnction to the 
calculation of von Neumann entropies and a definition of the orbital interac- 
tion is given in Section 2, which also contains a recipe for calculating these 
quantities using expectation values of operators in second quantisation. Sec- 
tion 3 presents the aspects of the DMRG relevant for this work. Subsequently, 
the orbital interaction is calculated for a set of test molecules in Section 4 in 
order to test the reliability of the proposed method. Section 5 then presents 
results for the ordering problem in the DMRG. Finally, the main findings are 
summarised in Section 6. 



2 Entanglement and orbital interaction 

2. 1 Definition of the orbital interaction Ip^g 

It is useful to divide the Hilbcrt space of a quantum-mechanical problem (the 
"universe") into two parts, which we will call the "system" and the "environ- 
ment" . If the basis of the system is described by a complete set of states {|i)} 
and that of the environment is described by the set the general wave 

function |^) can be written as [6] 

l^) = E^MK)b) • (1) 

With the help of the density operator for a pure state describing the total 
problem p = | \E' ) ( ^ | , one can express the reduced density matrix for the 
system as 

= E A,.'. = E 01 (^1 ^ )( * 1^') l^') = E cIjC,, . (2) 

This is the typical way to set up reduced density matrices. It requires the 
partition of the wave function as in Eq. (1) as well as the knowledge of the 
coefficients Cij. 

It is also possible, however, to formulate the problem in a different way. Sup- 
pose one has another, arbitrary basis for the system, then 

p:!:' = Epn,>', = E01(^l*)(*K)li) ■ (3) 

With the help of Eq. (1) one can write 

pZ' = T.C:,C,^,{^'\n'){n\^) = (*|P„>|*) , (4) 
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where 

^n> = Eli)K)H01 ■ (5) 

j 

Instead of a sum over the coefficients of the wave function, one can calculate 
the elements of the reduced density matrix with the help of expectation values 
of the operator Pn',n- The effect of Pn'^n is to change the state of the system 
from \n) to \n'). It can be viewed as a rotation in state space. 

After diagonalisation of p^^^ one obtains its eigenvalues uUa, which define the 
von Neumann entropy of the system [1,2], 

S^y^ ^ -Y^uj^lnu;^ . (6) 

a 

This quantity describes how much the system is entangled with the environ- 
ment for a given wave function |^). When the basis {|n)} describes only one 
orbital p, then the system contains only this orbital p and the corresponding 
entropy is the one-orbital entropy. If the system contains two orbitals p, q, 
then S^'^ is a two-orbital entropy. Since the two-orbital system is built up from 
two subsystems, namely, the orbitals p and q, one can apply the subadditivity 
property of S: 

SP1 <SP + S\ (7) 

where the equality holds when p and q are not entangled. The interpretation is 
straightforward: describes the entanglement of p and -S"* the entanglement 
of q with the rest of the environment, while S^"^ describes the entanglement of 
p and q with the rest of the environment. Any entanglement between p and q 
reduces S^"^ with respect to the sum of and S'^. Therefore, one can define 
the entanglement between two individual orbitals by 

I,,, = ^{S^ + S'^- S^'^) (1 - S,,) > , (8) 

where the Kronecker S ensures that Ip^p — 0, and the factor 1/2 prevents 
interactions from being counted twice. The quantity /p ^ is interpreted in the 
remaining part of this work as a measure of the orbital interaction. 

Although Ip^q describes the entanglement between two orbitals, it is not pos- 
sible to use Ip^q in order to build up entropies of larger systems. For example, 
when the system contains the orbitals p and the environment the orbitals q, 
then 

S'''> E ^P.. (9) 

pGsys,gGenv 

for all cases that we have investigated. The difference can amount to 60% of 

^sys 
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2.2 A recipe for the calculation S 



In order to calculate an orbital interaction using Eq. (8) , one has to determine 
and diagonalise the one- and two-orbital reduced density matrices and p^^. 
Their matrix elements are calculated in this work with the help of the right- 
hand side of Eq. (4), in which the matrix elements are written as expectation 
values of Pn',n with respect to the wave function I^P), where the operator 
changes the state of the system from \n) to \n'). 

One suitable way to represent the system basis {|^)} is the occupation-number 
representation. In particular, if the system consists of only one orbital p, 

{\n)} = {£j,T |0) '4a |0) ' |0) '4,t4 = {It) ' li) ' |0) , iTi)} (10) 

where c^ ^ are creation and annihilation operators for electrons with spin 
a in the orbital p, and |0) is the vacuum state. Eq. (10) leads to a one-orbital 
density matrix of dimension four. If the system contains two orbitals, then 
{\n)} consists of sixteen basis states, where the orbitals p and q are occupied 
by zero to four electrons 

{\n)} = {IT, 0) , IT, T) , IT, i) , IT, ti) , li, Ti) , • • • , ITi, Ti)} , (n) 

which leads to a 16 x 16 two-orbital density matrix. 

This representation determines not only the dimension of the resulting density 

matrices, but has two additional consequences. The first is that the matrix 
elements of the reduced density matrix pj^^^, are zero if |n) and \n') differ in 
the number of electrons or in the z-component of the spin, for example 

P|T).in> = E(*li)ITi)(TI01*) = o. (12) 

Since all of the basis states in Eq. (10) differ in either the number of electrons 
or the z-component of the spin, the one-orbital density matrix is diagonal in 
this representation, as shown in Table 1. The second consequence of Eqs. (10) 
and (11) is that the matrix elements (\E'| Pn',n |^) of Eq. (4) can be expressed 
using creation and annihilation operators. For example. 



pq 

^ITa>,|o,U> 



^E(^b')IT'i)(o,Til01*) 



(13) 



where hp^„ = c^p^^Cp^^ is the (occupation-) number operator for electrons with 
spin a in the orbital p. We exploit this property in order to determine the one- 
and two-orbital density matrices. 
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Table 1 

The one-orbital density matrix p^,f^, for an orbital p. The states {|n)}, {|n')} are 
defined in Eq. (10). The (occupation-)mimbcr operator for electrons with spin a in 

the orV)it ;til /; is , ' ' ' 



' P,(T 
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{np,i (1 - np,|)) 








(0| 








((1 - (1 - "p,t)) 





(Til 











("p,i^P,T) 



Consequently, the eigenvalues of in Table 1 are written in terms of the 
occupation-number operator. Thus, the resulting one-orbital entropy de- 
pends only on the occupancy of the orbital. In particular, = when the 
orbital is fully occupied or empty, and it reaches its maximal value = In 4 
when (fip^i) = (^p,t) = 0.5 and (^p,j.^p,t) = 0.25. 

In this work, the occupancies of the orbitals are the same for each spin (fip^i) — 
{np,-[), and, using the approximation {np^ifip,^) ~ (^p,j.) (^p,t)) S^ts 



S^np) 



(14) 



where Up = {hp^i) + {fip,]) is the occupancy of orbital p. Thus, S^lup) is a 
symmetric function of the occupancy in this case: SP{np) — 3^(2 — Up). We 
plot therefore the values of found in this work versus the deviation Aup 
of the occupancy from the Hartrec-Fock case in Fig. 1. One to six orbitals 
deviate from the approximation of Eq. (14), depending on the molecule. The 
deviations occur for Aup > 0.025. The largest one-orbital entropy found in 
this work < 0.31 is much smaller than the maximal value 1.38 
because the occupancies do not deviate strongly from the Hartree-Fock case 
and because {np,ihp,^) ^ {hp^i) {hp,f). 

Finally, in this representation it is possible to see that the two-orbital density 
matrices also contain the information contained in the reduced one-particle 
density matrix 
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and the reduced two-particle density matrix 

pS;° = E(^i^p>.,<^'|*) • 



(15) 



(16) 



The elements of p™^*^ are contained in the off-diagonal elements of the two- 
orbital density matrices (f^ while the elements of p*™° are contained in those 
diagonal elements of ff"^ where p and q are always occupied. 
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Fig. 1. One-orbital entropy as defined in Eq. (6) and Table 1 versus the deviation 
of the occupancy of orbital p, Up = {fip^i) + {fip^i), from the Hartree-Fock case 
{Aup = 2 — Up if Up > 1 and Aup = Up if Up < 1): circles - LiF, squares - CO, 
diamonds - N2, triangles - F2. Solid line - Eq. (14). The values of are taken from 
a DMRG calculation at m = 200 and Hartree-Fock ordering. 



This leads to the following recipe to calculate the orbital entanglement of 
Eq. (8): 

(1) Calculate an approximate wave function |^). 

(2) Calculate the one-orbital and two-orbital density matrices ff for every 
orbital p and q, and for every combination p, q. 

(3) Diagonalise the density matrices and calculate the entropies and 5*^^ 
using the eigenvalues uja- 

(4) Calculate Ip^q. 

In order to follow this recipe and to make use of the formulation of Eqs. (4) 
and (13), it is necessary to apply a method that can easily calculate expecta- 
tion values of operators in second quantisation with respect to the total wave 
function. This is the case for the DMRG, which is described in the following 
section. 



3 The DMRG and the calculation of von Neumann entropies 



This section presents only those aspects of the DMRG that are relevant for 
the following discussion. For further details we refer the reader to Refs. [3,7,8]. 
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3.1 The setup 



The DMRG is used here to determine the ground state and its energy for the 
time-independent, non-relativistic, electronic Hamiltonian (PuU-CI problem) 
with a well-defined error. The second-quantised Hamiltonian reads 

= XI '^P,q ^P,cr^q,o- + XI ^P,q,r,s ^l,cT^l,a''^r,a''^s,a ■ i^'^) 
PiSif p,q,r,s,a,a' 

Typically, T^^, Vp'^l^ are the one- and two-electron integrals in the basis of 
the canonical orbitals p, q, r, s, i.e., the eigenf unctions of the Fock-operator. 

The determination of all of the eigenstates of Eq. (17) requires a diagonalisa- 
tion in a Hilbert space which grows exponentially with the number of orbitals 
N: dim(7i) = 4^. The key idea of the DMRG is to form a reduced basis in 
an optimal, controlled way for the few eigenstates one is interested in, thus 
reducing the numerical effort. The size of the reduced basis B in this work is 
dim(B) = 16 • m^. 

For a given m, one determines the optimal basis iteratively with the help of 
the reduced density matrix (Eq. (2)). This requires the formation of a wave 
function as in Eq. (1) in every step. Consequently, one also must divide the 
total basis into a "system" and an "environment" part which is realised in the 
following way: First, one orders the orbitals onto a one-dimensional lattice, sec 
Fig. 2. Second, a boundary line defines a left and a right block. The product 
states that can be formed with the orbitals of each block are represented in 
this work by 4m states: the states of one block are the "system" states {|i)} 
and the states of the other block are the "environment" states {\j)} in Eqs. (1) 
and (2). 

In every step of the algorithm the boundary is moved by one site. When 

the separator is moved from left to right, the left block plays the role of 
the "system" and the right block the role of the "environment". When the 
separator is moved from right to left the roles arc exchanged. A zipper-like 
motion of the boundary back and forth through the whole lattice is called a 
sweep. 

Finally, it is important to point out that the ordering of the orbitals onto 
the lattice is in principle arbitrary. One can therefore choose to order the 
orbitals on the lattice with increasing energy, an ordering which we will call 
the Hartree-Fock (HF) ordering for the rest of this work. The infiuence on 
the energy convergence of different orderings, as mentioned in Section 1, is 
discussed in Section 5. 

In the subsequent sections, the ground states of LiF, CO, N2, and F2 are in- 
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Fig. 2. Setup of the DMRG: every orbital occupies one lattice site. In every step of 
the DMRG the separator is moved one site to the right or to the left. 

vestigated. The same cc-pVDZ basis set [9] has been used for all molecules, 
which gives rise to A?" = 28 canonical orbitals. We have calculated the one- 
and two-electron integrals using Dalton, a standard quantum-chemistry pro- 
gram package [10]. The molecules have been calculated at their experimen- 
tal distances: ruF = 156.3864 pm, rco = 112.8323 pm, rm = 109.768 pm, 
rF2 = 141.193 pm [11,12,13,14]. 



3.2 The error in energy and the measurement of Ip^, 



The states |i) of the system block are improved in every step of the algorithm. 
In particular, one diagonalises the Hamiltonian (Eq. (17)) in the reduced basis 
B, calculates the reduced density matrix p^^^ (Eq. (2)) and projects the 4m 
basis states representing the system block onto those m eigenstates of p^^/ that 
have the largest eigenvalues cua and are in that sense an optimal representation 
of the system block. The eigenvlues cUa also define the projection error of one 
step of the DMRG 

m 

PT"" = 1 - E (18) 

a=l 

since Y^a^a — 1- We take the projection error for a given m, P^, to be -P^^^ 
at the step when the separator is in the middle of the lattice in the last sweep, 
assuming that convergence in the number of sweeps has been achieved. 

If m is large enough so that dim^B = dimTY, then the FuU-CI problem is solved 
and Pm — 0. Since the DMRG is variational, every energy for a given m is 
larger than the exact energy E^, > -Eexact- Therefore, Em can be extrapolated 
to the exact (FuU-CI) ground state energy £'exact) which is in our case unknown, 
using [5] 

Em -E'cxact ~l~ (^P'rn ■ (1"^) 

In order to obtain inexact; we calculate a linear regression to Em{Pm) at six 
different m — (200,300,400,500,600). Typically, six sweeps are necessary to 
obtain a converged energy at m = 200 and four additional sweeps for each 
of the other values of m. Table 4 contains the extrapolated energies inexact- 
The errors given in Table 4 are the standard deviations for E^xact due to 
the extrapolation procedure, and we use error in this sense for the rest of this 
manuscript. There are other definitions in the literature for the error in energy 
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which also depend on P^^^ [15] and which consequently should give similar 
results. 

In order to calculate the orbital interaction 7p g, one has to form the oper- 
ators that lead to the reduced one- and two-orbital density matrices using 
Eqs. (4) and (11). For example, instead of evaluating Eq. (13) as a whole, one 
determines 



In this way, every matrix element of fl'^,^^, is expressed as sum of matrix el- 
ements of smaller operators. Only 23 of them are needed to construct (f^'^^t. 
They are stored in matrix form and must be transformed into the new basis at 
every step of the DMRG. We evaluate the 23 operators at the end of a DMRG 
calculation with the wave function \^) from the last step and calculate 
from the DMRG output. 

Although the one-orbital entropies can be determined throughout the DMRG 
algorithm by calculating the local densities {fii^a) via Table 1, it is necessary to 
also calculate them at the end of the procedure in order to obtain comparable 
accuracy in each matrix element and in order to calculate them with respect 
to the same wave function as the two-orbital entropies. 



4 The structure and accuracy of Jp ^ 

We have carried out calculations for the electronic ground state of the four test 
molecules (LiF, CO, N2, F2). We have obtained wave functions with different 
energies and errors by changing some parameters of the algorithm: the under- 
lying symmetry of the Hartree-Fock calculation, the ordering of the orbitals, 
and the parameter m. It turns out that the general structure of Jp^^ is not 
affected by these parameters. In other words, l^^q can already be determined 
using a comparatively inexpensive calculation, for example with m = 200 at 
HF ordering in any symmetry. 

The values of Ip^g are in the range < Ip^q < 0.18 in this work. The maximal 
values for the different molecules are Ip^q = 0.03 for LiF, Ip^g = 0.08 for CO, 
Ip^g = 0.11 for N2, and Ip^g = 0.18 for F2. The distribution of the Ip^g values is 
shown in Fig. (3). One can see that LiF, the only molecule in the series with 
an ionic bond, has a broad distribution of elements, while the other molecules 
show a somewhat smaller distribution and two or three intervals with a large 
concentration of Ip^g elements. The largest interaction Ip^g for F2 is one order 
of magnitude larger than the second largest interaction element. Also N2 and 
CO have two to three interaction elements that are clearly separated from the 
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Fig. 3. Histograms of l^^q calculated at HF ordering for m = 200: The number of Ip^q 
elements are counted for consecutive intervals. For example, the number of elements 
in 0.0001 < Ip,g < 0.0005 is displayed at 0.0005, the number in 0.0005 < ip,^ < 0.001 
is displayed at 0.001 and so on. 

rest, which is not the case for LiF. 

In order to visualise the structure of Jp^^, it is useful to assign a label to each 
orbital. In this work, the orbital labels stem from a Hartree-Fock calculation 
using the highest point group available in the Dalton program package: 
for F2, N2 and Civ for LiF, and CO. This means that the lowest- lying orbital 
in F2 always has the label la^, even for calculations in C\. Using these labels, 
we can specify the order in which the orbitals are put on the lattice in the 
DMRG calculation. The orderings are given in Tables 2 and 3, together with 
a label for the different cases. The labels designate the molecule, the ordering 
criterion, and the point group in which the HF calculation is carried out. For 
example, LiF-HF-Ci denotes a calculation for the LiF molecule, ordering by 
the orbital energies (HF ordering), and utilising a HF calculation in C\. 

The structure of 1^^^ can be examined using diagrams. We connect two orbital 
labels with a line if the corresponding value for 1^ ^ is larger than a chosen 
threshold, here l^^q > 0.005. In Fig. 4, we display such diagrams for the four 
molecules studied where the calculations have been carried out using the HF 
ordering and m = 200. From the diagrams, one can see the following features: 
First, in all cases the interaction couples predominantly orbitals of the same 
irreducible representation. Second, there are a few orbitals which are also cou- 
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Fig. 4. Diagram of Ip^q calculated at HF ordering and m = 200: Lines connect orbital 
labels with Ip^q > 0.01. The circle for CO and N2 denotes that the surrounding 
orbitals are all connected with each other. 

pled to orbitals of other irreducible representations (for example 3ai, 4ai, Ibi, 
lb2 in LiF). Third, those orbitals which couple different symmetry sectors are 
often energetically close to the highest occupied molecular orbital (HOMO). 
They have an occupancy Up that deviates from Up — 2 or Up — and thus 
have a large one-orbital entropy according to Section 2. 

Therefore, one can say that orbitals with a large one-orbital entropy also 
have large interactions Ip^q. They are often energetically close to the HOMO 
and are frontier orbitals in that sense (see Section 1). One can also ask, whether 
the value Ip g is connected to the integrals T^^ and Vp'^j.^ of the Hamiltonian 
in Eq. (17). To answer this question, we have a look at the largest Ip^g which 
always correspond to orbitals p — q which are az — cr*, tt^ — tt*, or tt^ — tt* 
combinations. The bonding and anti-bonding orbitals belong to the same irre- 
ducible representation for LiF and CO. For F2 and N2, however, these orbitals 
belong to different irreducible representations. Consequently, T^^ = for the 
latter cases and p, q are only connected by two-electron integrals in the Hamil- 
tonian, whereas for LiF and CO one can find one- and two-electron integrals. 
It is therefore not possible to deduce a clear correspondence between the size 
of Ip^q and the size of the one- and two-electron integrals. 

Another way to visualise the overall structure of Ip^g is to plot this quantity as 
a matrix for a given ordering. For example, the element 12,5 for case LiF-HF- 
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C2v denotes the interaction between the 2ai and Ibi orbitaL This also makes 
it possible to mirror the effect of different orderings: an ordering which groups 
strongly interacting orbitals together has large elements next to the diagonal 
in the plot of Ip^g. 

(a) (b) 
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Fig. 5. Ip^q calculated for ordering LiF-HF-C2j; (label defined in Table 2) with (a) 
m = 200 and (b) m = 600. 

In Fig. 5, we display plots of Ip^g for the case LiF-HF-C2^ with m — 200 and 
m — 600. The increase in accuracy has no significant effect on Ip g, although 
the electronic energies differ by about 7 • 10~^ a.u. (see Table 4). It can also 
be seen that the HF ordering results in large weight in the off-diagonal region 
of Ip^q and thus in an ordering for which strongly interacting orbitals are far 
apart. This holds also for all other cases of HF ordering. 

(a) (b) (c) 



q q q 




Fig. 6. Ip^g calculated for (a) N2-HF-£>2ft and m = 200, (b) N2-HF-C2^ and m = 600, 
and (c) N2-HF-C1 and m = 600 (labels defined in Table 3). 

In Fig. 6, one can see that the inffuence on Ip^q of different symmetries in the 
HF calculations is neghgible. The plots of Ip^q for N2 and HF ordering in D2/1, 
C2V: and Ci are almost identical. This is also refiected in Fig. 8 in which the 
curves of the electronic energies with respect to the DMRG steps lie on top of 
each other and consequently lead to the same extrapolated energy (see Table 
4)- 

Up to this point, the orbitals have been ordered on the DMRG lattice accord- 
ing to their energy (HF ordering). Different orderings result in different wave 
functions, energies, and structures of the Ip^g matrices. The next section deals 
with this issue in more detail. For now, only the effect of the different orderings 
on Ip^g is discussed. Fig. 7 shows plots of Ip^q for m — 600 and three cases from 
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(a) (b) (c) 

q q q 




5 1 0152 025 P 5 1 0152 025 P 5 1 0152 025 P 

Fig. 7. Ip^q from calculation using m = 600 and (a) F2-HF-D2/i, (b) F2-(23)-D2/j, 
and (c) F2-[18]-Z)2/i (labels defined in Table 3). The matrices in (b) and (c) are 
plotted in the HF ordering rather than the ordering of the actual calculation. 

Table 3: F2-HF-L'2/i, F2-(23)-D2/i, and F2-[18]-D2/i. There is a sizeable energy 
difference between the HF ordering and the other two cases of approximately 
6 • 10^^ a.u. (see Table 4). Despite these differences, no significant difference 
in the plots of Ip^q is discernible. In order to make this comparison possible, 
we have plotted all Ip^q matrices in the same ordering, namely the HF order- 
ing, although the actual calculations have been carried out using the labelled 
criteria. In Fig. 13 one can see plots of Ip^q for F2-(23)-D2/i and F2-[18]-D2/i in 
the ordering which has been used in the calculation. 

To conclude, we find that Jp ^ is a reliable definition of orbital interaction which 
can be calculated for small m in any ordering. The next section addresses the 
question of whether one can use the information in Ip^q to obtain an optimal 
ordering in the sense that the DMRG has a rapid convergence towards the 
exact wave function and energy. 



5 The ordering problem 

5.1 The problem and solution strategies 

Two things must be chosen before a DMRG calculation can be carried out: 
the basis and its ordering on the lattice. In principle, any orthonormal basis 
of orbitals can be used and the ordering is arbitrary. However, the ordering 
does affect the convergence of the DMRG in practice [4]. This can be seen by 
analyzing the behaviour of the energies as m is increased. With increasing m, 
the variational nature of the DMRG leads to a decrease in energy for every 
ordering. For a sufficiently large m, the difference to the exact energy can be 
made arbitrarily small. In this sense the ordering is arbitrary. However, the 
value of m that is necessary for a certain accuracy in energy depends on the 
ordering. For example, in Fig. 8 one can see that the case LiF-(23)-C2^ leads 
to a much lower energy for m — 200 than LiF-HF-C2t,. With increasing m 
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Table 2 

LiF, CO: Orderings of the orbitals as used in the DMRG: The labelling has 
the form (molecule) -(ordering criterion)-(point group) where (molecule) is LiF or 
CO, (ordering criterion) is HF (increasing orbital energy), (23) (using Eqs. (8), (21), 
(23)), or [18] (Ref. [18]). The label (point group) is the corresponding Schonflies 
symbol (for example C2v)- Orbital labels stem from calculations at the highest 
possible point group (here C2v)- Occupied orbitals are printed in bold face. 

Label Ordering 

LiF-HF-C2^ [lai 2ai 3ai 4ai Ibi Iba 5ai 2^2 26i 6ai 7ai 3bi 3^2 8ai 

9ai laa 46i 4^2 lOaillai 56i 5^2 12ail3ai 662 Gbi Uai 202] 

LiF-(23)-C2^ [2bi 36i 46i Ibi 56i 14ai 2a2 661 13ai 662 5ai 8ai 4ai lloi 

12ai3ail0ai 6ai lai la2 9ai 2ai 7ai 262 362 4^2 lb2 562] 

LiF-[18]-C2^ [lai2ai 9ai 7ai 5ai 6ai 14ail3ail0ail2ai 8ai 3ai lloi 4ai 

2bi 661 36i 46i 56i Ibi lb2 5^2 4^2 362 662 262 la2 102] 

C0-HF-C2^ [lai 2ai 3ai 4ai Ibi Iba 5ai 2bi 262 6ai 362 36i 7ai 8ai 

9ai 462 46i la2 lOai 56i 5^2 Ilail2ail3ai 2a2 662 661 14ai] 
CO-(23)-C2^ [56i 46i Ibi 2bi 5ai 7ai 36i la2 661 2ai Mailloi 9ai 4ai 

8ai 3ail2ai 6ai lai 2a2 13ai 662 5^2 4^2 lb2 262 3^2 lOoi] 

CO-[18]-C2^ [Iai2ail4ail3ail0ail2aillai 6ai 7ai 8ai 9ai 3ai 4ai 5ai 

662 562 362 462 262 lb2 Ibi 26i 46i 36i 5bi 661 la2 202] 

the energies for LiF-HF-C2,; improve and for m — 600 both cases show only a 
difference in energy of 3- 10~^ a.u., as seen in Table 4. However, the error in the 
energy in Table 4 is considerably larger for LiF-HF-C2,;. In the worst case, it is 
also possible to find orderings that cause a trapping of the DMRG algorithm 
in a local minimum [4]. The energy then does not decrease significantly with 
increasing m even though remains small. It is important to note that this 
trapping is not necessarily due to a certain ordering but is sometimes caused 
by an unfavourable choice of parameters for the DMRG calculation. 

In order to apply the DMRG to practical situations, it is therefore crucial 
to define a criterion for an optimal ordering for a given basis, in order to 
avoid trapping in local minima and in order to achieve the highest possible 
accuracy. The approaches which have been pursued so far have defined an 
orbital interaction and have ordered the orbitals so that strongly interacting 
orbitals are near each other on the one-dimensional lattice. 

One approach is to define the interaction between orbitals p and q in terms 
of the one- and two-electron integrals of the Hamiltonian [5,16]. Improved 
orderings then reduce the bandwidth of the T^^ matrix, for example. This 
approach is not able to avoid trapping in local minima [4] . In order to find an 
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Table 3 

N2 Orderings of the orbitals as used in the DMRG: The labeUing has the 
form (molecule)- (ordering criterion)-(point group) where (molecule) is N2 or F2, 
(ordering criterion) is HF (increasing orbital energy), (23) (using Eqs. (8), (21), 
(23)), or [18] (Rcf. [18]). The label (point group) is the corresponding Schonflies 
symbol (for example C2v)- Orbital labels stem from calculations at the highest 
possible point group (here D2h)- Occupied orbitals are printed in bold face. 

Label Ordering 

N2-HF-C1 [lag Ibiu 2ag 2biu Sag Ibgulbsu ^hg 1^2g 36i„ Aug 263^ 262n Sa^ 
~ C2v,D2h 2633 262p 46i„ Bbiu 6ag Ibig 363^ 362^ 66i„ la^ Tag 363^ 362^ 76i„] 
N2-(23)-i:>2fe[6ap Ibig 2b2g lb2g Ibsu 263^ 363^ la^ lag 76i„ 362u 5ag 2biu 3ag 

46iu 2ag 5biu 36i„ 362g Ibiu ^hg Tag 4.ag 2b2u lb2u ^hg "^hg 66i„] 
N2-[18]-i:)2/i[lag Tag 6ag 4.ag 5ag 2ag Sag la^ 363^ 263^ 163^ lb2g 2b2g Sb2g 

362u 2b2u Ib2ulb3u 263„ 363^ Ibig 2biu ^hu ^hu 56i„ 66iu Tbiu Ibiu] 
N2-[18]-C2t, [lag Ibiu Tbiu 66i„ Tag 5biu Gag Aag 36i„ bag 46i„ 2ag Sag 2biu 

3629 362« 263„ 263^ lb2g Ibaulbau 1^3^ 2623 262„ 363„ 363^ Ibig la„] 
N2-[18]-Ci [Ibiu 3629 66i„ Tag Ibig 362« 262« 36i„ 263^ bag 2ag 2biu 1633 lb2u 

Ibsu lb2g Sag Abiu 2b2g 263^ Aag 363^ 6ag 56iu la^ 76iu 363^ lag 
F2-HF-D2^ [lag Ibiu 2ag 2biulb2ulb3u 3ag Ib2glb3g36i„ 263^ 2b2u 46i„ 2625 

263^ Aag Sag 5biu 6ag 362u 363^ Tag Ibig 66i„ la„ 362g 363^ 76i„] 
F2-(23)-i:>2/i[5ap 2biu 4ag Sag 36i„ 46i„ 6ag Ibiu lag 363^ 362« 263^ Ibsg 262^ 

lb2u 2ag bbiu 76i„ 66i„ 7ag 16ig la„ 3623 363„ Ibsu 263^ lb2g 2629] 
F2-[18]-L'2h [lag Tag 6ag bag Aag 2ag Sag 36i„ 2biu 46iu bbiu 66i„ 76i„ Ibiu 

3629 2629 Ib2glb3g 2633 3639 ^hu 2fe3u Ib3ulb2u 252u 362u I&I9 lOu] 

ordering criterion based on or V^j^lj, a recent study [17] has investigated 
a large number of orderings for the Cr2 molecule using a genetic algorithm, 
which has not yet led to a general criterion for different molecules and basis 
sets. 

Another approach is to group the orbitals according to their irreducible repre- 
sentations and then order the orbitals within these groups in order to maximise 
the one-orbital entropy along the lattice [4,18]. The net effect is that some 
entangled (interacting) orbitals are placed close together but are also some- 
what distributed over the lattice. This is called "competition between entan- 
glement localisation and interaction localisation" in Ref. [18]. The label "[18]" 
is attributed to this criterion, which obviously depends on the underlying sym- 
metry of the HF calculation. This is illustrated for N2-[18]-£'2/i, N2-[18]-C2^, 
and N2-[18]-Ci in Fig. 9. The distribution of along the lattice is plotted for 
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Table 4 

Energies and errors due to different orderings: Label - (molecule)- (ordering 
criterion)- (point group) defined in Tables 2 and 3; m - size of the largest used Hilbert 
space dim(;B) = 16 ■ m^; energies - electronic energy in atomic units (without nuclei 
interaction), extrapolation according to Eq. (19) for m = 200, 300, 400, 500, 600 error 
is the standard deviation. 
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the three point groups investigated. The maxima correspond to boundaries 
between groups of orbitals with the same irreducible representation. For the 
Ci ordering in Fig. 9(a), the sites with large one-orbital entropy are bunched 
up in the middle of the lattice. In Fig. 8, we display the effect of these order- 
ings on the convergence of the DMRG: only in the case N2-[18]-C2t, can one 
see an improvement over the HF ordering N2-HF-Z)2/j- Therefore, one is not 
guaranteed that a new ordering according to the strategy of Ref. [18] leads to 
improved energy convergence. 
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DMRG sweeps 



Fig. 8. Electronic energy (no nucleus-nucleus interaction) versus number of DMRG 

sweeps: m = 200 (sweeps 1 to 6), m = 300 (sweeps 7 to 10), m = 400 (sweeps 11 to 
14), m = 500 (sweeps 15 to 18), m = 600 (sweeps 19 to 22): solid line - HF ordering; 
open symbols - ordering using Eqs. (8), (21), (23); filled symbols - ordering using 
Ref. [18]; circles - C2v', squares - D2h; triangles - Ci 



(a) (b) (c) 




Lattice Site 

Fig. 9. Orbital entropies SP following Eq. (6) and Table 1 for (a) N2-[18]-Ci, (b) 
N2-[18]-C2^, and (c) N2-[18]-D2h (labels defined in Table 3). 

5.2 New strategy using Ip^q 

A strategy based on the orbital interaction Ip^g combines the advantages of 
the two earlier approaches. On the one hand, one is using the entanglement 
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information of a many-body wave function as in Ref. [18]. On the other hand, 
one has a specific interaction in matrix form as in Refs. [5,16,17]. The quan- 
tity is not sufficiently specific because it only indicates how much orbital 
p interacts with all the other orbitals, while Ip^q is a direct measure of the 
interaction between p and q. The analysis of Ip^q in Fig. 4 has also shown that 
orbitals from the same irreducible representation are coupled strongly. It is 
therefore not necessary to account for this fact separately as had to be done 



Here we search for an improved ordering which localises the interaction Ip^q, 
i.e., reduces the bandwidth of the Ip g matrix. The optimal ordering is found 
using simulated annealing [19]. We argue that this approach is an improvement 
over the Cuthill-McKee algorithm, used for example in Refs. [5,16], which dis- 
tinguishes only between occupied and unoccupied matrix elements and thus 
neglects the information contained in the size of the elements. The annealing 
algorithm is comparable to the genetic algorithm [17] because both probe dif- 
ferent, randomly generated configurations. The cost function in the annealing 
process plays a similar role to the fitness function in the genetic algorithm. 

Following the idea that good orderings should place strongly interacting or- 
bitals near each other, we have investigated several cost functions F that all 
favour orderings in which large elements of Ip^q are on the secondary diagonal 
and which have a small bandwidth. It turns out that a cost function 



where r = |j9 — g|, is a good starting point. In order to increase the attraction 
to the secondary diagonal of large elements of Ip^q, we set 



We have also made further adjustments to the cost function F. From the 
discussion of the ordering criterion in Ref. [18], we have learned that strongly 
interacting orbitals should not be bunched up, (for example, in the middle 
of the lattice) but should instead be more evenly distributed. Therefore, we 
set r = 0.5 for elements on the secondary diagonal in regions of length A^/5 
around the edges and the middle of the lattice, where N is the number of 
orbitals, i.e.. 



in Ref. [18]. 




0.5 if |p - g| = 1 
|p — g| otherwise . 



(22) 



r = < 



0.5 if \p — q\ — 1 and 

{{p, q} < N/5 or {p, q} > N - N/5 or 
N/2 - N/10 < {p, q} < N/2 + N/10) 

|p — g| otherwise . 



(23) 
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An ordering created by the use of Eqs. (8), (21) and (23) will be labelled 
"(23)". 



In all cases in which orderings have been generated using Eq. (23) we have 
found an improved energy convergence with respect to calculations with HF 
ordering, and the calculations are not trapped. However, we find that there 
is no clear correspondence between the value of the cost function for a given 
ordering and the respective energy convergence. In other words, small changes 
in the ordering and the cost function can influence the energy convergence 
severely. For example, the case LiF-[18]-iI'2/i has a similar energy convergence 
as the case LiF-(23)-D2fe as shown in Table 4. However, F = 0.6397 in the 
former case and F = 1.574 in the latter, while the HF ordering yields F = 
0.281. 



To summarize, while, in our opinion, Eqs. (21) and (23) represent an improved 
cost function, one cannot always be certain that an ordering determined by 
these equations really is the optimal one. While the value of F contains valu- 
able information, additional information, which can be obtained by the visual 
investigation of Ip^g as a matrix, is also important. This will be discussed for 
speciflc cases in the following. 
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Fig. 10. Ip^q calculated at m 
defined in Table 2). 



600 for (a) LiF-(23)-C2^, and (b) LiF-[18]-C2^ (labels 



LiF, for example, is a rather unproblematic case: any new ordering that re- 
duces the bandwidth of Ip^g results in improved convergence in the DMRG. In 
Fig. 10(a), one can see that the bandwidth of Ip^g for the ordering of this work 
is reduced compared to the HF ordering, which results in a reduced extrap- 
olated energy and a reduced error (see Table 4 and Fig. 8). While the result 
using the criterion of Ref. [18] looks more compact (Fig. 10(b)), it yields simi- 
lar results for the energy convergence. A comparison of the orbital orderings in 
Table 2 shows that both cases group together orbitals of the same irreducible 
representation. 

For CO and N2, the situation is more complicated. In Fig. 11(a), one can see 
that the form of Jp g for our ordering is again more spread out than for the 
ordering of Ref. [18] (Fig. 11(b)). This time, however, the latter criterion leads 
to a more exact energy, as can be seen in Fig. 8. 
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Fig. 11. /p,y calculated at m = 600 for (a) CO(23)-C2^, and (b) CO-[18]-C2^ (labels 
defined in Table 2). 
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Fig. 12. calculated at m = 600 for (a) N2-(23)-L>2h, (b) N2-[18]-D2h, (c) 
N2-[18]-C2^', and (d) N2-[18]-Ci (labels defined in Table 3). 

Fig. 12(a) shows Jp.q for N2 plotted for the ordering criterion of this work. 
Figs. 12(b)-(d) are determined by the criterion of Ref. [18] and use different 
symmetries in the underlying HF calculation. All four plots have a reduced 
bandwidth compared to the HF ordering shown in Fig. 6, but only the cases in 
Figs. 12(a) and 12(c) show better energy convergence. The criterion of Ref. [18] 
yields a slightly more compact form of /p g, as can be seen in Fig. 12(c), and 
leads to a better convergence in energy than our criterion (see Fig. 8). The 
cases displayed in Figs. 12(b) and 12(d) distribute or accumulate the inter- 
acting orbitals too much, leading to poor energy convergence, as explained in 
the discussion of Fig. 9. This underscores that subtle changes in the structure 
of the Ip^q matrix influence the energy convergence of the DMRG. 

For F2, in Fig. 13 one cannot say that one plot of Ip^q is more compact than 
the other. Despite this similarity, the differences in the energy convergence are 
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Fig. 13. calculated at m = 600 for (a) F2-(23)-r>2ft and (b) F2-[18]-r>2ft (labels 
defined in Table 3). 

as pronounced as for CO. Only this time the criterion of our work leads to the 
more exact energy. 

To conclude, we can say that the apphcation of Ip^q, Eqs. (21) and (23) leads to 
an ordering with a considerably better energy convergence than a HF ordering, 
and the results do not depend on the symmetry used for the underlying HF 
calculation. In addition, orderings with a good energy convergence have the 
following properties: the bandwidth of the Ip^q matrix is small, large elements 
/p q are grouped on the secondary diagonal, and pronounced accumulation and 
scattering of large Ip^q elements are avoided. 

However, we have not been able to establish a distinct correspondence between 
the orbital interaction and the energy convergence, although Ip^q is a very reli- 
able quantity. Therefore, one should test a few different orderings determined 
by varying some of the parameters of the annealing process before one sets 
up a DMRG calculation aimed at high accuracy. The resulting energy conver- 
gence can be checked for small sizes of the Hilbert space. This is demonstrated 
in Fig. 8, where the orderings with the lowest energy for m = 200 also have 
the lowest energy for m — 600. 

5.3 Change of the basis 

Evidently, it is desirable to obtain additional insight into the mechanism for 
energetic convergence in the DMRG. One way to achieve this is not only 
to consider the ordering of the orbitals but also the choice of the orbitals 
themselves. Since the canonical orbitals are not a mandatory choice, it is 
possible to construct a new basis that might suit the DMRG better. So far, 
there has been one attempt to use a localised basis in the DMRG [20] that has 
not, however, led to an improved convergence relative to canonical orbitals. 

An obvious choice are natural orbitals, i.e., the eigenfunctions of the one- 
particle density matrix (Eq. (15)). They lead to rapid convergence of the 
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configuration interaction scheme, and should therefore be favourable for the 
DMRG as well. Since the one-particle density matrix is contained in p^'', one 
can easily construct approximate natural orbitals, for example, from an ap- 
proximate wave function of a calculation at m = 200 with HF ordering. 



The natural orbitals must also be ordered on the lattice. One can use an order- 
ing according to occupation number as a reference ordering. DMRG results for 
the electronic energies for the four sample molecules calculated with canonical 
and natural orbitals for various orderings are displayed in Fig. 14. One can 
see that there is a small increase in convergence for the natural orbitals com- 
pared to the canonical orbitals in the HF ordering in every case. One can then 
apply the ordering criterion of this work which consistently improves energy 
convergence. However, the optimal energy convergence does not exceed that 
found for the best orderings of canonical orbitals (not shown in Fig. 14). 
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Fig. 14. Electronic energies (no nucleus-nucleus interaction) versus number of 
DMRG sweeps for m = 200 (sweeps 1 to 6). Solid line: canonical orbitals and HF 
ordering, dashed line: natural orbitals ordering by eigenvalues, dashed line (circles, 
squares): natural orbitals ordering with Eqs. (8), (21), (23) (6*2^, D2h)- 

Finally, one can consider an iterative improvement of the natural orbitals: 
Beginning with a DMRG calculation using canonical orbitals in HF ordering, 
one can obtain an improved wave function which can then be used to calculate 
natural orbitals which in turn can serve as a basis for a next DMRG calculation 
yielding a wave function with an even lower energy. This procedure can be 
repeated yielding an improved wave function and energy at every iteration 
until the energy converges. We have applied this iterative procedure and have 
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found that the energies decrease only for two or three iterations and then start 
to fluctuate. Therefore, further investigation is needed to understand how to 
construct an optimal basis for the DMRG. 



6 Summary 

In this work, we have used concepts from quantum information theory to 
formulate a definition of orbital interaction Ip^q. For a given wave function, 
Ip^q is defined by the subtraction of the entanglement of two orbitals taken 
together with the rest of the system from the sum of the entanglement of 
two individual orbitals with the rest of the system. The advantage of this 
definition is that one includes information beyond the Hartree-Fock treatment. 
The disadvantage is that a correlated wave function must be calculated. Given 
a correlated wave function, we have developed a recipe in Section 2.2 for the 
calculation of Ip^g. This recipe additionally provides an alternative method to 
calculate one-orbital entropies which are also of central importance in other 
work [4,18]. 

We have calculated Ip^q using correlated wave functions obtained from a DMRG 
calculation and the recipe of Section 2.2. The resulting interaction Ip g does 
not depend strongly on the accuracy of the underlying wave function. For the 
four test molecules we have treated, we have shown that it is possible to cal- 
culate Ip^q in a reliable fashion. The structure of Ip^q is consistent with another 
criterion based on one-orbital entropies and with chemical intuition: partially 
occupied orbitals of the same irreducible representation interact strongly. 

As an application, we have used Ip^q to study the ordering problem in the 
DMRG, in which one has to order orbitals on a one-dimensional lattice so 
that strongly interacting orbitals are near each other. We have developed a cost 
function for a simulated annealing process that leads to an improved ordering 
for all of the cases we have treated, i.e., the subsequent DMRG calculation 
leads to a lower energy. We have also found that orderings with a good energy 
convergence have a small bandwidth in the Ip^q matrix, and have a distribution 
of large matrix elements on the secondary diagonal that is relatively uniform 
along the lattice. However, we have not been able to identify a consistent 
scheme that leads to an optimal ordering of the orbitals for the DMRG with 
this approach. Careful checks and additional attempts to find better orderings 
are needed. 

A more general solution to this problem might lie in the construction of the 
basis itself. We have therefore investigated the influence of the use of natural 
orbitals on the convergence of the DMRG. This leads to a slight improvement 
over the use of canonical orbitals, but the ordering problem still remains. It 
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is known, however, that the DMRG yields excellent results with a basis of 
Pz orbitals for conjugated polymers. Hence, an optimal basis should consist 
of orbitals that are localised in real space and are close in energy. It might 
then be possible to construct an optimal basis for the DMRG for which the 
ordering is either obvious or irrelevant. Then the ordering problem should be 
of minor importance. This will be the topic of a subsequent study. 
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